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1  .Introduction 


Equilibrium  problems  for  many  physical  systems  are  modelled  by 
parameter  dependent  nonlinear  equations 

(1.1)  F(z,X)  =  0.  1  .  . 

ch.  -  ' 

Under  fairly  general  conditions  the  set  of  solutions  (z,V)  of  fMi  forms  a 
differentiable  manifold,  and,  typically  in  applications,  interest  centers  not  so 
much  on  computing  a  few  solutions  but  rather  on  analysing  the  form  and  special 
features  of  this  manifold.  Mn  this  paper  w^-rctefflify  some  of  the  sources  of  the 
errors  which  are  necessarily  arising  in  such  a  computational  analysis.-.. 


For  the  presentation  we  adopt  the  same  setting  as  in  [18]  or  [33].  Let  X  ^ 
and  Y  be  real  Hilbert  spaces  and  F:  X  m  Y  a  Fredholm  mapping  of  class 
r>2,  and  index  p  £  1  for  which  the  domain  contains  an  open  set  S  of  X.  A  point 
x  e  X  is  regular  if  the  Jacobian  DF(x)  maps  X  onto  Y.  Then  it  is  well  known  that 
the  set  of  all  regular  solutions, 


(1.2) 


M  =  { x  ;  x  e  S,  F(x)  =  0,  x  regular),  , 


./  / 


r  , 

is  a  p-dimensional  Cf-manifold  in  X  without/  boundary.  By  restricting 
consideration  to  this  regular  solution  manifold  we  have  assumed  that  a  suitable 
unfolding  of  the  problem  has  been  chosen,(see  e.g.  [20]).  The  tangent  space 
TXM  at  xe  M  may  be  identified  with  the  keiftelof  DF(x);  that  is, 


(1.3) 


TXM  ■  ker  DF(x)  =  [u;  u  e  X,  DF(x)u  =  0 }, 


1  This  work  was  in  part  supported  by  the  Office  of  Naval  Research  under 
contract  N-00014-80-C-9455  and  the  National  Science  Foundation  under  grant 
DCR-8309926. 
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whence 

(1.4)  NXM  =  (TxM)1  =  rgeDF(x)* 

is  the  normal  space  at  the  point1.  A  given  p-dimensional  subspace  T  of  X 
induces  a  local  coordinate  system  of  M  at  x  e  M  if 

(1.5)  To  NXM  ={0}. 

More  specifically,  if  (1.5)  holds  then  there  exist  open  neighborhoods  Vi  and  V2 
of  the  origins  of  T  and  X,  respectively,  as  well  as  a  unique  CM  -  function 
wiVi-^T1  ,w(0)  =  0,  such  that 

(1.6)  M  n  V2  =  {  y  e  X  ;  y  =  x  + 1  +  w(t),  t  e  V^, 

(see  e.g.  [33]).  If  (1.5)  holds  then  x  is  a  non-singular  point  with  respect  to  the 
given  coordinate  space  T,  else  it  is  a  singular  point  or  foldpoint  with  respect  to 
T.  Clearly,  any  point  xe  Mis  regular  with  respect  to  its  the  tangent  space  TXM  . 

As  indicated  already  by  the  form  of  (1.1),  many  applications  involve  a 
natural  orthogonal  splitting 

(1.7)  X  =  Z©A,  Z  =  A1- 

of  the  domain  space  X  into  a  state  space  Z  and  a  p-dimensional  parameter 
subspace  A.  Then  interest  centers  on  determining  the  singular  points  with 
respect  to  A  .  In  equilibrium  problems  these  may  be  expected  to  be  the  points 
where  a  change  in  the  stability  behaviour  of  the  physical  system  occurs. 

The  basic  procedures  for  the  computational  analysis  of  such  a  manifold 
M  are  the  continuation  methods.  When  M  has  dimension  p  >  1 ,  these  methods 
require  a  restriction  to  some  path  on  M  and  then  produce  a  sequence  of  points 
along  that  path.  Obviously,  it  is  not  easy  to  develop  a  good  picture  of  a  multi¬ 
dimensional  manifold  solely  from  information  along  such  paths.  This  led 


1  As  usual,  the  asterisk  denotes  the  Hilbert  space  adjoint. 


recently  to  the  development  of  methods  for  the  computation  of  simplicial 
approximations  of  p-dimensional  subsets  of  M.  Two  intrinsically  different 
methods  of  this  type  have  been  presented  in  [1],  [2]  and  [35],  [36]. 

Besides  these  methods  for  computing  certain  sets  of  points  on  M,  another 
important  class  of  numerical  procedures  concerns  the  detection  and 
computation  of  foldpoints  on  M  with  respect  to  a  given  coordinate  space  T,  such 
as  the  natural  parameter  space  A.  There  is  a  large  and  growing  literature  about 
such  methods;  further  references  can  be  found  in  [23],[24],[26]. 

This  presentation  will  address  only  some  aspects  of  these  two  classes  of 
procedures.  There  are,  of  course,  further  related  computational  tasks.  For 
example,  we  may  wish  to  compute  contour  lines  or  contour  surfaces  of  a 
prescribed  functional  on  M.  Alternately,  many  applications  lead  to  the  need  for 
computing  certain  solutions  of  some  differential  equation  on  a  manifold.  This  is 
equivalent  with  the  solution  of  a  differential-algebraic  system  of  equations 
which,  in  the  case  of  ordinary  differential  equations,  is  a  topic  with  its  own 
burgeoning  literature. 

Our  interest  here  is  not  to  elaborate  on  these  various  computational 
techniques  but  in  identifying  some  of  the  errors  arising  in  their  connection.  More 
specifically,  in  the  next  section  we  discuss  the  error  induced  by  a  discretization 
of  the  basic  mapping  F.  Then,  in  Section  3  we  turn  to  error  questions  connected 
with  continuation  methods  while  the  final  Section  4  identifies  related  questions 
in  the  computation  of  simplicial  approximations  of  parts  of  M  as  well  as  of 
foldpoints  on  the  manifold. 


2.  Discretization  Errors 

In  practice,  the  mapping  F  often  represents  a  differential  operator 
involving  several  parameters  which  then,  for  the  computation,  has  to  be 
approximated  by  a  finite-dimensional  analogue.  More  specifically,  suppose  that 
we  have  a  natural  splitting  (1.7)  of  the  domain  space  X  and  that,  as  in  (1.1),  our 
problem  is  written  in  the  form 


(2.1) 


F(z,X)  =0 


ze  Z  ,  *K  €  A. 


Then  only  the  state  variable  z  needs  to  be  discretized;  in  other  words,  the 
approximating  equations  will  have  the  form 

(2.2)  Fh(Zh,X)=0,  zeZh,  U  A, 

where  now  Fh  maps  a  discretized  space  Xh  =  Zh  ©  A  to  another  such  space  Yh. 
Since  the  parameter  space  is  unchanged  we  may  expect  that  the  regular 
solutions  of  the  discretized  problems  form  a  p-dimensional  manifold  Mh  in  Xh- 

Frequently,  Xh  can  be  imbedded  into  X  which  then  allows  for  a  direct 
comparison  of  the  two  manifolds  M  and  Mh.  However,  in  order  to  measure  the 
distance  between  these  two  manifolds,  and  hence  the  discretization  error,  we 
need  to  choose  some  coordinate  system.  More  specifically,  suppose  that  at  the 
desired  point  x  e  M  a  local  coordinate  system  has  been  induced  by  the  p- 
dimensional  subspace  T  of  X.  In  other  words,  in  a  neighborhood  of  x  the  points 
of  M  are  represented  as  x(t)  =x+t+w(t),  te  T,  w(t)eT1,  w(0)=0.  If  the 
approximation  is  sufficiently  close,  then  some  part  of  Mh  should  belong  to  the 
domain  of  validity  of  this  local  coordinate  system  and  hence  there  should  be  a 
unique  point  xh  on  Mh  which  has  the  same  coordinate  t=0  as  x;  that  is,  which  is 
in  x+T1 .  We  may  expect  also  that  T  induces  a  local  coordinate  system  at  Xh  on 
Mh,  and  hence  that,  locally  near  Xh,  the  points  of  Mhare  representable  in  the 
form  Xh(t)  =Xh+t  +Wh(t)  ,  t  eT,  Wh^eT1,  Wh(0)  =  0.  Now  it  makes  sense  to 
measure  the  discretization  error  as  the  distance  ||x(t)-Xh(t)||  in  X  between 
points  on  the  two  manifolds  which  have  the  same  local  coordinate  t  e  T.  In  other 
words,  the  discretization  error  is  a  strictly  local  concept  and  depends  on  the 
choice  of  the  particular  local  coordinate  system. 

If  a  natural  parameter  decomposition  (1.7)  is  available  and  the 
coordinate  system  is  induced  by  the  choice  T  =  A  ,  then  we  compare  points  with 
the  same  X-values  and  the  discretization  error  measures  the  difference 
between  the  states  of  these  points.  This  is  often  proposed  as  the  definition  of 
the  discretization  error.  But,  clearly,  it  is  only  a  feasible  choice  as  long  as  A 
induces  a  local  coordinate  system  at  the  point  x  e  M  ;  that  is,  as  long  as  x  is  not 
foldpoint  with  respect  to  A.  Since  such  foldpoints  are  of  central  importance  in 
many  applications,  this  is  certainly  not  a  generally  acceptable  approach. 
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The  development  of  a  rigorous  theory  of  these  discretization  errors  is  a 
fairly  recent  undertaking.  For  one-dimensional  solution  manifolds  a  priori 
estimates  were  first  developed  in  [12]  and  then  [14].  The  latter  results  were 
generalized  in  [16]  to  manifolds  of  arbitrary  dimension.  In  the  three  parts  of  [12] 
the  cases  of  non-singular  points,  limit  points,  and  simple  bifurcation  points  were 
treated  separately  and  specific  local  coordinates  were  used  at  each  of  these 
points.  In  [14],  [16]  general  local  coordinate  systems  in  the  above  sense  are 
utilized  to  avoid  the  requirement  of  distinguishing  between  the  different  types  of 
points. 


All  these  results  involve  a  family  of  approximate  problems  (2.2)  which 
converge  in  some  sense  to  the  original  problem  (2.1)  when  the  real 
discretization  index  h  >  0  tends  to  zero.  More  specifically,  the  results  in  [14], [16] 
are  based  on  projection  methods;  that  is,  a  family  {Ph}  of  finite-rank  projections 
is  assumed  to  be  given  on  the  range  space  Y  for  which  limh_>oPh  y  =  y  for  each 
y  e  Y.  Then,  with  a  bounded  linear  operator  A  from  X  to  Y  such  that  ker  A  =  A, 
and  Az=A|Z  is  an  isomorphism  onto  Y,  we  can  define  the  approximate 
mappings 


Fh  :Xh  — »  Yh,  Xh  =  Zh©A,  Yf,  =  PhY,  Zh  =  Qh  Z  ,  Qh=  A^PhAz. 

Now,  if  a  particular  stability  condition  holds,  the  estimate 

Co  II  (l-Ph)Ax(t)  ||  ^  ||  x(t)  -  xh(t)  ||  £  Ci  ||  (l-Ph)Ax(t)  ||  , 

is  valid  for  all  sufficiently  small  t  e  T  and  h  >  0  with  constants  Ci>Co>0  that  are 
independent  of  h  and  t.. 

A  different  approach  was  taken  in  [17].  There,  only  a  single  discretized 
equation  is  considered  instead  of  a  converging  family  of  such  equations,  and 
estimates  are  obtained  which  correspond  to  the  local  error  estimates  in  the 
numerical  solution  of  initial  value  problems  for  ordinary  differential  equations. 
These  estimates  have  been  applied  to  an  analysis  of  the  behavior  of  the 
socalled  reduced  basis  techniques  pioneered  in  structural  mechanics  (see 
[15], [28]). 
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These  a  priori  estimates  are  of  considerable  theoretical  interest,  but  for 
the  practical  application  we  require  a  posteriori  estimates  which  measure  the 
error  of  the  specific  computed  points  on  the  approximate  solution  manifold  Mh. 
Such  estimates  are  needed  not  only  to  judge  the  reliability  of  the  computed 
results  but  also  to  control  adaptive  procedures  aimed  at  achieving  dependable 
results  within  a  given  range  of  accuracy  at  minimal  cost. 

For  finite  element  discretizations  there  exists  a  growing  literature  on  the 
computation  of  such  a  posteriori  estimates;  we  refer  only  to  [4]  and  the  two 
proceedings  [3],  [6]  where  many  other  references  can  be  found.  These  results 
concern  largely  linear  problems,  but  more  recently  they  are  also  being 
extended  to  nonlinear  problems,  see  e.g.,  [5],  [8], [32],  [34],  For  example,  the 
approach  in  [34]  supposes  that  --  in  the  notation  of  the  a  priori  estimates  --  the 
local  coordinate  space  T  belongs  to  the  approximate  domain  space  Xh.  If  the 
desired  exact  point  x  and  the  approximate  computed  point  Xh  are  sufficiently 
close,  then  it  can  be  shown  that 

(2.3)  ||x-  xh||  =  ( 1  +o(1)  )||y||  as  ||x  -  xh||  ->  0 


where  y  is  the  solution  of  the  linearized  equation 
(2.4)  F(xh)  +  DwF(xh)y  =  0,  W  =  T1 . 


The  discretization  of  (2.4)  has  the  exact  solution  yh  =  0  and  hence  ||y||  is  the 
discretization  error  of  (2.4)  which  can  be  estimated  by  means  of  one  of  the 
known  a  posteriori  estimators  for  linear  equations. 


Various  computational  examples  for  such  a  posteriori  estimates  in  the 
case  of  finite-element  discretizations  of  one-dimensional  manifolds  have  been 
given  in  the  cited  articles.  Corresponding  estimates  for  two-dimensional 
manifolds  are  used  in  the  nonlinear  adaptive  finite  element  solver  NFEARS 
developed  jointly  by  I.Babuska,  C.K.Mesztenyi,  and  W.C.Rheinboldt.  These 
results  indicate  that  the  effectivity  of  these  error  estimates  corresponds  to  that  for 
linear  problems  and  that  the  computational  cost  is  only  a  fraction  of  the  cost  of 
l  solving  the  nonlinear  problem.  In  general,  the  errors  vary  considerably  along  a 

!  solution  manifold  which  shows  the  importance  of  adaptive  mesh-refinements. 

» 
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However,  there  remain  many  open  problems.  In  fact,  it  appears  to  be 
reasonable  to  require  that  such  estimates  (a)  apply  to  large  classes  of 
problems,  discretizations  and  error  norms,  (b)  are  dependable  and  cost- 
effective  in  the  range  of  engineering  accuracies,  (c)  allow  for  the  use  of  most 
applicable  local  coordinate  system,  and  (d)  provide  also  estimates  of  the  error 
of  computed  foldpoints.  At  this  time  we  are  still  very  far  from  such  type  of 
estimates,  and  accordingly  it  is  well  justified  to  characterize  the  computation  of  a 
posteriori  estimates  for  nonlinear  parametrized  problems  as  being  in  an  early 
stage  of  development. 

As  stressed  before,  all  these  various  error  estimates  are  local  in  nature 
and  hence  cannot  provide  information  about  the  comparative  global  shapes  of 
the  exact  and  approximate  manifolds.  In  particular,  it  turns  out  that  the  two 
manifold  may  have  a  different  number  of  connected  components.  In  other 
words,  there  may  be  numerical  solutions  which  do  not  approximate  exact 
solutions.  Such  spurious  solutions  have  been  observed  in  many  contexts  (see 
e.g.  [9], [19], [27]). 

A  very  simple  example  for  such  spurious  solutions  arises  in  connection 
with  the  boundary  value  problem 

(2.5)  u"  +Xsinu  =  0  ,  0  <  s  <  1,  u(0)  =  u(1)  =  0 

which  can  be  viewed  as  a  model  of  the  classical  Euler  rod.  For  X  <  0,  there 
exists  only  the  trivial  solution  u  =  0,  which  is  physically  plausible  since  under 
pure  tension  the  rod  should  remain  straight.  Nevertheless,  discretizations  of 
(2.4)  constructed  either  by  finite  differences  or  by  finite  element  techniques 
exhibit  non-zero  solutions  for  certain  negative  values  of  X  (see  e.g.  [5]  for  some 
picture).  These  particular  spurious  solutions  are  not  especially  disturbing  since 
they  occur  in  a  region  of  little  interest  and  disappear  to  infinity  when  the  mesh  is 
refined. 

Unfortunately,  not  all  spurious  solutions  are  so  "harmless".  For  instance, 
the  following  problem 


(2.5) 


u"  +  a  (1-u)  exp(-Xy(1+u))  =  0  ,0S  t  £1,  u(0)  =  u(1)  =  0 


models  an  exothermic  chemical  reaction  in  a  slab.  It  can  be  shown  that  all 
solutions  must  have  the  symmetric  property  u(t)  =  u(1-t),  t  e  [0,1].  However,  a 
discretization,  such  as 

(2.6)  -Xj-i+2xj  -Xj+i  =  h  a  (1  -xj)  exp(-X/(1+Xj))  ,  i=1,...,n-1 ,  x0  =  xn  =  0,  h  =  1/n, 

models  the  same  reaction  only  inside  each  one  of  a  collection  of  n  cells  while 
across  the  interfaces  there  is  merely  mass  transport.  As  a  consequence 
unsymmetric  solutions  of  (2.6)  exist  which  branch  off  from  the  (discretized) 
symmetric  solution  at  certain  symmetry  breaking  bifurcation  points.  These 
branches  tend  to  the  exact  symmetric  solution  when  h  goes  to  zero,  but,  clearly, 
they  do  not  correspond  to  any  similar  feature  of  the  original  problem,  (see  e.g. 
[10],  [11]). 

For  practical  applications  it  is  certainly  of  considerable  interest  to  provide 
error  estimators  which  identify  spurious  solutions.  In  [5]  it  was  shown  that,  for  a 
problem  such  as  (2.4),  the  spurious  solutions  carry  error  estimators  which  are 
much  larger  than,  say,  any  tolerance  acceptable  in  engineering  calculations 
and  hence  which  can  be  rejected  on  that  basis.  There  is  considerable  need  for 
a  more  detailed  study  of  this  question. 

3.  Continuation  Methods 

As  noted  earlier,  the  basic  procedures  for  the  computational  analysis  of 
our  manifold  M  are  the  continuation  methods  and  these  methods  require  an  a 
priori  restriction  to  some  path  on  M.  In  other  words,  if  M  has  dimension  p  >  1 , 
then  our  (discretized)  equations  must  be  augmented  by  p-1  suitable  equations 
which  specify  the  desired  path  on  M.  As  a  consequence  the  continuation 
methods  are  always  applied  to  an  equation  of  the  form 

(3.1 )  F0(x)  =  0  ,  F0:  Rn  -»  R"-i  ,  n£2 

involving  an  operator  F0  which  is  of  class  Cr,  r  >2,  on  some  open  set  in  Rn.  We 
denote  its  regular  solution  manifold  by  N  =  {x  e  Rn;  Fo(x)=0  ;  rankDF0(x)=n*1}.  It 
is  readily  verified  that,  on  the  set  of  regular  points  of  F0,  the  mapping 
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(3.2)  u:  Rn— ►  R".  DF0(x)u(x)  =  0,  ||u(x)[|2  =  1 .  det  (DF0(x)Tu(x))  >  0 

defines  a  CM-vector  field.  Let  x°  e  N  be  a  given  point  and  y:J  ->  N  the  unique, 
maximally  extended  solution  of  the  initial  value  problem 

(3.2)  y'(s)  =  u(y(s))  ,  V  s  e  J  ,  y(0)  =  x°, 

then  the  trajectory  y(J)  is  exactly  the  connected  component  N0  of  N  containing 
x°  (see,  e.g.,  [33]).  This  connection  between  the  initial  value  problem  (3.2)  and 
its  first  integral  (3.1)  has  been  used  by  numerous  authors. 

All  continuation  processes  begin  from  a  given  point  x°  e  M  and  produce 
a  sequence  of  points  xk,  k=0,1 ,2,...,  which  approximate  points  of  N.  For  any  k>0, 
the  step  from  xk  to  xk+1  corresponds  to  an  implementation  of  the  local 
coordinate  representation  (1 .6).  More  specifically,  if  T  =  span(t),  t  e  Rn,  t  *  0,  is 
a  local  coordinate  space  of  N  at  xk  ,  then,  for  any  fixed  y  e  Rn,  the  Jacobian  of 
the  augmented  equations 

(33)  (tVy))  ■  0 

is  non-singular  for  all  x  in  some  neighborhood  of  xk.  Thus,  if  y  approximates  a 
point  of  N  in  that  neighborhood,  then  it  follows  readily  that  (3.3)  has  a  unique 
solution  xk+1  e  N  which  can  be  computed  by  means  of  a  locally  convergent 
iterative  process  applied  to  (3.3)  and  started,  say,  at  y. 

The  various  continuation  methods  differ  in  the  choice  (a)  of  the 
coordinate  space  T,  (b)  of  the  predicted  point  y,  and  (c)  of  the  local  iterative 
process  for  solving  (3.3).  In  most  cases  a  linear  prediction  y  =  xk+hv  is  chosen, 
whence  (b)  subdivides  into  the  choice  (bl)  of  the  predictor-direction  v  e  Rn,  and 
(b2)  of  the  step-length  h  >  0.  The  so-called  pseudo-arclength  methods  use  for  t 
and  v  the  tangent  vector  u(xk)  of  N  at  xk  ,  while  in  the  PITCON  code  (see  [30]), 
only  the  prediction  v  =  u(xk)  is  along  the  tangent  and  T  is  specified  by  a 
suitable  natural  basis  vector  t  =  e',  i  =  ik,  of  Rn.  On  the  other  hand,  in  the 
PLTMG  code  [7]  both  t  and  v  are  a  linear  combination  of  u(xk)  and  a  specified 
vector  characteristic  for  the  problem.  The  local  iterative  process  (c);  that  is.  the 
corrector,  usually  is  a  chord-Newton  method  with  the  Jacobian  of  the  mapping 


(3.4)  at  xk  or  y  as  the  iteration  matrix.  Other  corrector  methods  include,  in 
particular,  the  multigrid  approaches  and  some  chord  Gauss-Newton  method. 

In  the  construction  of  the  predicted  point  it  is  highly  desirable  to  estimate 
the  interval  of  validity  of  the  local  coordinate  system  at  xk.  This  is  equivalent  with 
the  requirement  of  estimating  the  distance  from  xk  to  the  nearest  foldpoint  of  N 
with  respect  to  the  chosen  local  coordinate  system.  At  present,  there  are  no 
reliable  methods  for  that  purpose.  A  second  critical  question  concerns  the  error 
||y-xk+1||  of  the  prediction  and  many  processes  indeed  utilize  estimates  of  this 
error.  For  example,  in  PITCON  such  an  estimate  is  obtained  as  the  difference 
between  the  chosen  linear  predictor  and  some  quadratic  extrapolation  involving 
in  turn  an  approximation  of  the  local  curvature.  This  is  then  applied  in  the 
construction  of  the  step  length  h. 

Whenever,  the  corrector  applied  to  (3.4)  converges  to  a  point  where  the 
Jacobian  remains  non-singular,  then  this  limit  must  be  a  point  of  N  But  no 
estimate  of  the  prediction  error  by  itself  can  guarantee  that  the  next  point  xk+1 
will  again  belong  to  N0  and  not  to  some  other  component  of  N.  An  important 
control,  which  can  be  employed  here,  is  the  orientation  of  N0  given  by  (3.2). 
This  is  indeed  incorporated  in  PITCON  and  has  proved  capable  of  identifying 
many  of  the  potential  "jumps"  of  the  process  between  components  of  N.  But,  of 
course,  there  remain  cases  where  jumps  to  components  with  identical 
orientation  are  possible.  Clearly  there  is  still  much  room  for  other  estimates  of 
the  prediction  error  and  for  further  controls  that  signal  when  the  process  has  left 
the  connected  component  N0. 

Since  all  solutions  of  (3.4)  with  non-singular  Jacobian  are  points  of  N 
(although  not  necessarily  of  N0),  the  quality  of  the  approximation  of  N  by  the 
computed  points  {xk}  is  controlled  by  the  size  of  the  termination  error  of  the 
corrector.  The  estimation  of  this  termination  error  has  been  the  topic  of  many 
studies,  but  from  the  convergence  theory  of  iterative  processes  for  nonlinear 
systems  of  equations  it  is  well-known  that  exact  bounds  for  these  errors  require 
information  about  the  mapping  in  a  sufficiently  large  neighborhood  of  the 
solution  (see  e.g.  [13]) .  Moreover,  these  bounds  do  not  take  into  account  the  ill- 
conditioning  introduced  by  round-off.  In  fact,  it  is  well-known  that  the  uncertainty 
region  of  the  solution  of  a  system  as  (3.4)  can  be  quite  large.  For  example,  if  the 
quartic  polynomial 


x4  -  202  x3  +  1 529  x2  -  51 4898  x  +  6497400  =  0,  x  e  Ft1 , 


with  the  roots  x  =  49,50,51 ,52  is  evaluated  by  means  of  the  standard  Horner 
scheme  in  rounded  8-decimal  digit  arithmetic,  then  the  size  of  the  uncertainty 
interval  {x;  |x-52|  <  5}  of  the  largest  root  is  about  5=  0.09  ,  (see,  e.g.,  [31]). 
Clearly,  the  occurence  of  ill-conditioning  of  such  a  magnitude  may  destroy 
completely  the  reliability  of  the  overall  process.  The  control  of  this  phenomenon 
appears  to  be  one  of  the  most  promising  areas  for  the  application  of  interval 
techniques  in  connection  with  continuation  methods. 

4.  Triangulations  and  Foldpoint  Calculations 

As  observed  already  in  the  Introduction,  continuation  processes  require 
us  to  develop  a  picture  of  the  p-dimensional  manifold  (1.2)  solely  from 
computed  information  along  a  priori  specified  paths.  The  method  introduced  in 
[35],  [36]  for  obtaining  simplicial  approximations  (triangulations)  of  open  subsets 
of  the  manifold  uses  almost  the  same  tools  as  these  continuation  methods  and 
hence  has  also  a  similar  error  behavior. 

Suppose  that  our  problem  is  defined  by  the  finite  dimensional  mapping 
(4.1)  F:  Rn  ->  Rm,  n=  m+p,  p>2 

and  that  M  denotes  its  regular  solution  manifold.  The  triangulation  process  then 
begins  with  the  choice  of  a  simplicial  decomposition  Z  of  RP,  such  as,  for 
instance,  the  well-known  Kuhn-triangulation,  or,  for  p=2,  a  triangulation  of  R2  by 
means  of  equilateral  triangles.  The  aim  is  to  transfer  the  knots  of  some  part  of  Z, 
together  with  their  connectivity  information,  from  RP  onto  M.  As  in  any 
continuation  methods,  a  starting  point  x  on  the  manifold  M  is  assumed  to  be 
known.  In  its  basic  form  the  process  then  consists  of  two  steps:  First  a  suitable 
"patch"  of  the  reference  triangulation  Z  is  mapped  onto  the  affine  tangent  space 
x+TxM.  using  an  appropriate  basis  of  TXM.  Thereafter,  a  locally  convergent 
iterative  process  is  applied  to  "project"  the  resulting  knots  of  the  mapped 
simplices  from  x+TxM  onto  the  manifold  M.  These  two  steps  are  then  repeated 
with  one  of  the  computed  points  on  M  in  place  of  the  original  point  x.  But,  of 


course,  knots  of  I  that  have  already  been  mapped  onto  M  will  not  be  used 
again. 
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In  order  to  ensure  that  the  images  of  the  simplices  of  I  on  M  form  a 
simplicial  approximation  of  a  portion  of  the  manifold,  we  need  to  choose  bases 
on  the  tangent  spaces  TXM  which  change  smoothly  from  point  to  point.  In  other 
words,  we  require  a  moving  frame  on  some  open  subset  M0  of  the  manifold  M. 
Recall  that  a  moving  frame  of  class  Cs,  0<s<r,  associates  with  any  x  €  M0  an 

ordered  basis  (frame)  (u1 . uP}  of  TXM  such  that  each  coordinate  map  u',  1<i<p, 

from  M0  into  the  tangent  bundle  TM  defines  a  vector  field  of  class  Cs  on  M0. 
Hence,  we  need  an  algorithm  which  generates  for  each  point  x  of  M  an  n  x  p 
matrix  U(x)  with  orthonormal  columns  such  that  DF(x)U(x)  =  0  and  that  the 
mapping  U:  M0  RPxn  is  of  class  Cs  on  M0.  Standard  methods  for  computing 
tangent  bases  generally  do  not  produce  continuously  varying  matrices  U(x). 

A  moving  frame  algorithm  of  the  desired  form  has  been  introduced  in 
[35],  It  uses  an  n  x  p  reference  matrix  Tr  with  orthonormal  columns  and  assumes 
that  some  method  is  available  for  computing  at  any  xe  Man  orthonormal  basis 
matrix  U(x)  of  TXM  which,  of  course,  is  not  expected  to  depend  continuously  on 
x.  Let  M0  be  the  open  subset  of  M  where  the  subspace  T  of  Rn  spanned  by  the 
columns  of  Tr  induces  a  local  coordinate  system.  If  x  e  M0  and  we  compute  the 
singular  value  decomposition  A(x)T(U(x)Tr)B(x)  =  D(x),  then  it  was  shown  in  [36] 
that  the  mapping  x  -+  U(x)A(x)B(x)T  is  of  class  CM  on  M0  and  hence  defines 
the  desired  moving  frame.  If  the  dimension  p  of  the  manifold  is  small  in 
comparison  with  the  space  dimension  n,  then  the  principal  cost  of  this  algorithm 
derives  from  the  computation  the  original  basis  matrix  U(x). 

There  are  various  ways  of  implementing  the  computation  of  the  tangent 
bases  U(x)  and  of  the  local  iterative  process.  For  instance,  in  the  cited  articles 
the  QR-factorization 

DF(x)T  -  0(p)  . 

gives  U(x)  as  the  last  p  columns  of  the  orthogonal  n  x  n  matrix  Q,  and  the  chord 
Gauss-Newton  process,  defined  by 


(i)  solve  RTz  *  F(yk)  forzeRP,  (ii)  yk+1  *  yk-Q(z,0)T,  k=0,1 ,2 


converges  locally  to  the  unique  point  y*  e  M  with  the  coordinate  t  =  y°-x  in  the 
local  coordinate  system  induced  by  TXM.  On  the  other  hand,  in  NFEARS  the 
original  equations  are  augmented,  as  in  PITCON,  by  appropriate  unit  basis 
vectors  of  Rn,  and  the  tangent  basis  is  derived  from  a  triangular  factorization. 
Then  a  chord  Newton  method  is  applied  as  the  local  iterative  process. 

There  are  two  controlling  aspects  in  this  process.  One  of  them  concerns 
the  behavior  of  the  local  iterative  process  and  hence  raises  analogous 
problems  as  for  continuation  processes.  As  before,  three  basic  issues  are 
again  the  estimation  of  the  prediction  error  and  of  the  termination  error,  as  well 
as  the  control  of  the  effects  of  round-off.  In  addition,  the  earlier  problem  about 
the  interval  of  validity  of  the  local  cocordinate  system  becomes  here  a  question 
about  the  numerical  behavior  of  the  moving  frame  algorithm.  In  particular, 
whenever  we  apply  the  algorithm  at  one  of  the  already  computed  knots  x  on  M, 
we  should  be  able  to  determine  whether  this  point  x  still  belongs  to  the  open 
subset  M0  of  M  where  the  subspace  T  of  Rn  spanned  by  the  reference  matrix  Tr 
induces  a  local  coordinate  system.  This  leads  again  to  the  question  of 
estimating  the  distance  from  x  to  the  set  of  foldpoints  of  M  with  respect  to  T.  As 
mentioned  earlier  there  exists  no  method  for  this.  Fortunately,  numerical 
evidence  suggests  that,  in  connection  with  the  triangulation  process,  this  lack  is 
often  not  very  critical.  A  related  open  problem  concerns  the  quality  of  the 
computed  moving  frame  when  we  come  close  to  the  boundary  of  M0. 

As  noted  in  the  Introduction,  in  many  application  interest  centers  on 
computing  the  foldpoints  of  M  with  respect  to  a  natural  parameter  space  A.  But 
space  limitations  allow  us  only  to  give  here  a  very  cursory  summary  of  some  of 
the  basic  issues  as  they  relate  to  our  present  setting.  As  defined  earlier,  a  point 
x  e  M  is  a  foldpoint  with  respect  to  A  if  the  intersection  An=AoNxM  is  non-trivial; 
that  is,  equivalently,  if  there  is  a  non-zero  tangent  vector  of  x  which  belongs  to 
the  state  space.  The  integer  q  =  dim(Afsi)  is  the  first  singularity  index  of  x.  With 
the  natural  splitting  (1.7)  we  can  write  our  equation  in  the  form  (2.1).  Then  a 
foldpoint  x  =  (z.A.)  e  M  has  first  singularity  index  q  exactly  if  the  range  in  Y  of  the 
partial  derivative  DzF(x)  of  F  has  co-dimension  q.  The  cut  Mo(x+An1)  of  M 
through  x  orthogonal  to  An  characterizes  the  type  of  the  foldpoint  (see  e.g.  [18]). 

We  mention  here  only  a  few  of  the  numerous  computational  problems 
associated  with  foldpoints: 


(a)  Detection  problem:  Given  a  set  {x1 . xk}  of  points  in  X  and  appropriate 

additional  information  about  the  problem  at  these  points,  determine  whether 
some  foldpoint  x*  of  M  is  in  a  neighborhood  of  this  set. 

(b)  Approximation  error:  For  a  given  approximation  x  e  X  of  a  foldpoint  x*  of  M, 
estimate  the  distance  ||  x  -  x*||  in  X  between  these  two  points. 

(c)  Foldpoint  computation:  For  foldpoints  x*  of  M  of  a  prescribed  type  design 
iterative  processes  which  are  locally  convergent  to  x*  from  any  sufficiently  close 
approximation  x  e  X. 

(d)  Discretization  error:  If  Xh  is  a  computed  foldpoint  of  the  discretized  equations 
(2.2)  which  approximates  a  foldpoint  x  of  the  original  problem  (2.1),  compute  an 
a  posteriori  estimates  of  the  distance  ||x  -  Xh|| . 

(e)  Foldset  triangulation:  Compute  a  simplicial  approximation  of  a  set  of 
foldpoints  of  M  of  a  prescribed  type. 

Not  only  the  problems(b)  and  (d),  but  also  all  the  other  problems  involve 
error  estimation  questions.  For  example,  problem  (a)  requires  in  essence  an 
inclusion  procedure  for  the  zeros  of  certain  nonlinear  mappings.  This  is  most 
easily  seen  in  the  case  of  a  mapping  (3.1)  underlying  a  continuation  method 
where  we  assume  that  the  natural  splitting  (1.7)  is  given.  Then,  theoretically,  the 
foldpoints  of  the  regular  solution  manifold  N  are  the  points  x  e  Rn  which  solve 
the  constrained  problem 

<4-2>  (det(D°zFo(x)) )  -  0  •  rank(DF0(x))  -  n-1 . 

Evidently,  on  the  basis  of  information  at  finitely  many  points  alone  we  cannot 
expect  to  obtain  a  guaranteed  inclusion  result  for  (4.2);  at  least  some 
information  about  F0  in  a  suitable  neighborhood  of  the  given  set  of  points  is 
needed  here.  This  is  already  seen  for  very  simple  inclusion  results.  Suppose 
that  A  =  span{v}  and  that  x1  and  x2  are  two  successive  points  produced  by  a 
continuation  procedure.  If  both  points  belong  to  the  same  connected  component 
of  N,  and  ui  =  u(xi),  j=1,2,  are  the  corresponding  oriented  tangent  vectors  (3.2), 
then  sgn^u1)  *sgn(vTu2)  implies  that  there  is  a  limit  point  of  odd  order 
between  the  two  points.  Evidently,  without  further  information  about  F0  we 
cannot  verify  that  both  points  belong  to  the  same  connected  component  of  N 


and  without  that  assumption  the  statement  need  not  hold.  On  the  other  hand, 
suppose  that  we  only  know  that  x1  ,x2  e  N  and  that  we  orient  the  tangent  vectors 
by  means  of  (u1)T(u2)  >  0.  Now,  if  the  determinants  det(DF0(xi)T,ui),  j=1 ,2,  have 
different  signs  and  the  points  are  sufficiently  close,  then  they  cannot  belong  to 
the  same  connected  component  of  N,  and  hence  there  must  be  a  bifurcation 
point  of  odd  order  "between"  them.  These  observations  can  be  extended  readily 
to  the  more  general  mappings  (4.1)  using  points  xi  e  M,  j=1,...,p+1,  which,  for 
instance,  form  the  corners  of  a  p-simplex.  But,  even  without  entering  into  any 
details  this  discussion  already  indicates  that  procedures  for  the  general 
problem  (a)  indeed  require  information  about  the  behavior  of  the  mapping  on 
appropriate  sets  and,  in  general,  also  a  restriction  to  specific  types  of  foldpoints. 

Problem  (b)  depends  on  the  measure  of  the  distance  between  the 
foldpoint  x*  and  its  approximation.  Suppose  that  we  consider  again  the  finite 
dimensional  case  (4.1);  then  the  size  of  the  smallest  principal  angle  between 
the  m-dimensional  normal  space  rge(DF(x))  and  the  p-dimensional  natural 
parameter  space  A  gives  some  information  about  the  approximation.  These 
principal  angles  are  readily  computed,  but  again  further  information  about  the 
mapping  is  needed  to  derive  from  this,  say,  an  estimate  ||x-x*||.  Except  for 
simple  special  cases,  there  appears  to  exists,  at  present,  no  reliable, 
computable  estimator  for  this  norm  difference. 

As  noted  in  the  introduction,  there  is  a  large  literature  on  local  iterative 
methods  for  computing  certain  types  of  foldpoints.  Many  of  these  methods 
begin  with  a  suitable  augmentation  of  the  given  equation  such  that  the  Jacobian 
of  the  resulting  system  is  square  and  non-singular.  This  allows  then  the 
application  of  a  more  or  less  standard  iterative  process  for  their  solution. 
Obviously,  it  would  lead  to  far  to  enter  here  into  any  details  of  the  numerous 
augmentations  that  have  been  considered.  But,  clearly,  in  all  these  cases  the 
behavior  of  the  process  is  once  again  essentially  controlled  by  the  error  of  the 
initial  approximation,  the  termination  error  of  the  method,  and  the  influence  of 
round-off  on  the  reliability  of  the  result.  It  should  also  be  noted  that  for  specific 
types  of  foldpoints  there  are  other  approaches  which  utilize  the  particular 
characterization  of  these  points.  In  particular,  for  limit  points  on  one-dimensional 
manifolds,  as  they  arise  in  connection  with  continuation  methods,  there  exists  a 
wide  variety  of  such  procedures,  as,  for  instance,  the  comparative  study  [25] 
shows. 
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As  noted  in  Section  2,  a  computed  foldpoint  of  the  discretized  problem 
(2.2)  need  not  approximate  a  corresponding  foldpoint  of  the  original  problem 
(2.1).  But  if  it  does,  then  it  is  certainly  of  considerable  interest  to  estimate  the 
distance  between  these  two  points.  There  appear  to  exist  only  very  few  results 
along  that  line.  We  mention  here  merely  the  a  posteriori  estimates  given  in  [5] 
for  the  case  of  simple  limit  points  and  for  a  finite  element  discretization  of  a  two- 
point  boundary  value  problem.  Further  results  in  the  important  area  of  problem 
(d)  are  certainly  very  much  needed. 

Finally,  in  connection  with  problem  (e)  we  mention  only  the  methods  in 
[22]  and  [29]  for  computing  sequences  of  points  along  certain  paths  of  limit 
points  on  two-dimensional  regular  solution  manifolds.  All  these  methods  are 
closely  related  to  the  continuation  procedures  and  hence  involve  essentially 
the  same  error  questions.  There  appear  to  exist,  at  present,  no  methods  for 
computing  simplicial  triangulations  of  higher  dimensional  sub-manifolds  of 
foldpoints  of  any  type. 
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